Simulating Tetrathionate Circuit with Bioscrape and BioCRNpyler

Liana Merk

BE240 Spring 2020

In [1]:
from bioscrape.types import Model
from bioscrape.sbmlutil import import_sbml
from bioscrape.simulator import py_simulate_model

from biocrnpyler import *

import numpy as np
import pylab as plt
import tqdm
import pandas as pd

import matplotlib as mpl
import holoviews as hv
import bokeh.plotting
import bokeh.io
from bokeh.models import Range1d
from bokeh.layouts import gridplot
bokeh.io.output_notebook()
hv.extension('bokeh')

import seaborn as sns
sns.set()
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/__init__.py:23: UserWarning: No module named 'dnaplotlib'
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/__init__.py:24: UserWarning: plotting is disabled because you are missing some libraries
Loading BokehJS ...

BioCRNpyler Custom Classes

Phosphorylation

In [2]:
class Phosphorylation(Mechanism):
    def __init__(self,name="phosphorylation",mechanism_type="phosphorylation"):
        Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
    def update_species(self,s1, **keywords):
        phosphoprotein = ComplexSpecies([s1,Species("P",material_type="phosphate")])
        return [phosphoprotein,Species("P",material_type="phosphate")]

class AutoPhosphorylation(Phosphorylation):
    def __init__(self,name="autophosphorylation",mechanism_type="phosphorylation"):
        Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
    def update_reactions(self,s1,component = None, kphos = None, kdephos = None, \
                                                        part_id = None,**keywords):
                                                        
        if part_id == None:
            part_id = s1.name+"-P"
        if kphos == None and component != None:
            try:
                kphos = component.get_parameter("kphos",part_id = part_id,mechanism = self)
            except ValueError:
                kphos = 0
        if kdephos == None and component != None:
            try:
                kdephos = component.get_parameter("kdephos",part_id = part_id,mechanism = self)
            except ValueError:
                kdephos = 0
        #if component == None and (kphos == None or kdephos == None):
        #    raise ValueError("Must pass in a Component or values for kphos, kdephos.")
        if((kphos == 0 and kdephos == 0) or (kphos == None and kdephos == None)):
            print(component)
            print("kphos is {} and kdephos is {}".format(kphos,kdephos))
            raise ValueError("Must pass in at least one value, kphos or kdephos for "+part_id)
        rxns = []
        phosphate = Species("P",material_type="phosphate")
        phosphoprotein = self.update_species(s1)[0]
        if(kphos == None):
            kphos = 0
        if(kdephos == None):
            kdephos = 0
        if(kdephos > 0 and kphos > 0):
            rxns += [Reaction([s1,phosphate],[phosphoprotein],k = kphos,k_rev=kdephos)]
        elif(kphos > 0):
            rxns += [Reaction([s1,phosphate],[phosphoprotein],k = kphos)]
        elif(kdephos > 0):
            rxns += [Reaction([phosphoprotein],[s1,phosphate],k = kdephos)]
        
        return rxns
class SimpleKinaseMechanism(Phosphorylation):
    def __init__(self,name="Kinase",mechanism_type="phosphorylation"):
        Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
        self.phosphate = Species("P",material_type="phosphate")
    def update_species(self,s1,s2):
        #s1 causes s2 to get phosphorylated. So, we need to update species
        #with s2's phosphorylated form
        return [ComplexSpecies([s2,self.phosphate]),self.phosphate]
    def update_reactions(self,s1,s2,component = None, k = None,  \
                                                    part_id = None,**keywords):
        if part_id == None:
            part_id = s1.name+"-"+s2.name
        if k == None and component != None:
            k = component.get_parameter("kkinase",part_id = part_id,mechanism = self)
        s2_p = ComplexSpecies([s2,self.phosphate])
        return [Reaction([s1,s2],[s1,s2_p],k=k)]

class TwoStepKinaseMechanism(Phosphorylation):
    def __init__(self,name="Kinase",mechanism_type="phosphorylation"):
        Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
        self.phosphate = Species("P",material_type="phosphate")
    def update_species(self,s1,s2):
        #s1 causes s2 to get phosphorylated. So, we need to update species
        #with s2's phosphorylated form
        s2_p = ComplexSpecies(s2,self.phosphate)
        s1_s2_complex = ComplexSpecies([s1,s2])
        s1_s2_plus_p = ComplexSpecies([s1,s2_p])
        return [s2_p,s1_s2_complex,s1_s2_plus_p,self.phosphate]
    def update_reactions(self,s1,s2,component = None, kb1 = None, ku1=None,kkinase=None, kb2 = None, ku2=None, \
                                                    part_id = None,**keywords):
        """this mechanism involves three reactions:
        1) s1 binds to s2 to make Complex(s1,s2) this uses kb1 and ku1 for binding and unbinding
        2) Complex(s1,s2) + phosphate makes Complex(s1,Complex(s2,phosphate)) this uses kkinase, and only goes forwards
        3) Complex(s1,Complex(s2,phosphate)) unbinds to form s1, Complex(s2,phosphate). this uses kb2 and ku2
        """
        if part_id == None:
            part_id = s1.name+"-P-"+s2.name
        if kkinase == None and component != None:
            kkinase = component.get_parameter("kkinase",part_id = part_id,mechanism = self)
        if ku1 == None and component != None:
            ku1 = component.get_parameter("ku1",part_id = part_id,mechanism = self)
        if kb1 == None and component != None:
            kb1 = component.get_parameter("kb1",part_id = part_id,mechanism = self)
        if kb2 == None and component != None:
            kb2 = component.get_parameter("kb2",part_id = part_id,mechanism = self)
        if ku2 == None and component != None:
            ku2 = component.get_parameter("ku2",part_id = part_id,mechanism = self)
        specs = self.update_species(s1,s2)
        binding = Reaction([s1,s2],specs[1],k=kb1,k_rev=ku1)
        kinase = Reaction([specs[1],self.phosphate],[specs[2]],k=kkinase)
        kunbinding = Reaction([specs[2]],[s1,specs[0]],k=ku2,k_rev=kb2)
        return [binding,kinase,kunbinding]

class PhosphoTransfer(Phosphorylation):
    def __init__(self,name="phosphotransfer",mechanism_type="phosphorylation"):
        Mechanism.__init__(self,name=name, mechanism_type=mechanism_type)
    def update_reactions(self,s1,s2,component = None, k = None,  \
                                                        part_id = None,**keywords):
        if part_id == None:
            part_id = s1.name+"-"+s2.name
        if k == None and component != None:
            k = component.get_parameter("ktransfer",part_id = part_id,mechanism = self)
        s1_p = self.update_species(s1)[0]
        s2_p = self.update_species(s2)[0]

        return [Reaction([s1_p,s2],[s2_p,s1],k=k)]

Combinatorial Binding

In [3]:
class Combinatorial_Cooperative_Binding(Mechanism):
    """a reaction where some number of binders bind combinatorially to a bindee"""
    def __init__(self,name="Combinatorial_Cooperative_binding",
                            mechanism_type="cooperative_binding"):
        Mechanism.__init__(self,name,mechanism_type)
    def make_cooperative_complex(self,combo,bindee,cooperativity):
        """given a list of binders and their cooperativities, make a complex
        that contains the binders present N number of times where N is
        each one's cooperativity"""
        complexed_species_list = []
        for binder in combo:
            binder_cooperativity = int(cooperativity[binder.name])
            #I hope that cooperativity is an int! what if it isn't
            if(binder_cooperativity > 1):
                complexed_species_list += [Multimer(binder,binder_cooperativity)]
            else:
                complexed_species_list += [binder]
        complexed_species_list += [bindee]
        if(len(complexed_species_list)==1):
            myspecies = complexed_species_list[0]
        else:
            myspecies = ComplexSpecies(complexed_species_list)
        return myspecies
    def update_species(self,binders,bindee,cooperativity=None,\
                            component = None, part_id = None, **kwords):
        cooperativity_dict = {}
        for binder in binders:
            binder_partid = part_id+"_"+binder.name
            if ((cooperativity == None) or (type(cooperativity)==dict and binder_partid not in cooperativity) \
                                                                                        and (component != None)):
                coop_val = component.get_parameter("cooperativity", part_id = binder_partid, mechanism = self)
            elif type(cooperativity)==dict and binder_partid in cooperativity:
                coop_val = cooperativity[binder_partid]
            if component == None and ( cooperativity == None):
                raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
            cooperativity_dict[binder.name]=coop_val
        #allbinders = binders+[bindee]
        out_species = []
        for i in range(1, len(binders)+1):
            for combo in it.combinations(binders,i):
                #go through every possible combination of reactants and dna and make
                #all the complexes
                out_species += [self.make_cooperative_complex(combo,bindee,cooperativity_dict)]
        return out_species
    
    def update_reactions(self,binders,bindee,component=None,kbs=None,kus=None,\
                                                            part_id = None,cooperativity=None,**kwords):
        binder_params = {}
        for binder in binders:
            binder_partid = part_id+"_"+binder.name
            if ((type(kbs)==dict and binder not in kbs) or (type(kbs)!=dict and component != None)):
                kb = component.get_parameter("kb", part_id = binder_partid, mechanism = self)
            elif(type(kbs)==dict and binder in kbs):
                kb = kbs[binder.name]
            elif(type(kbs)!=dict and component == None):
                raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
            if ((type(kus)==dict and binder not in kus) or (kus == None and component != None)):
                ku = component.get_parameter("ku", part_id = binder_partid, mechanism = self)
            elif(type(kus)==dict and binder in kus):
                ku = kus[binder.name]
            elif(type(kus)!=dict and component == None):
                raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
            if ((cooperativity == None) or (type(cooperativity)==dict and binder.name not in cooperativity)  \
                                                                                        and component != None):
                coop_val = component.get_parameter("cooperativity", part_id = binder_partid, mechanism = self)
            elif type(cooperativity)==dict and binder.name in cooperativity:
                coop_val = cooperativity[binder.name]
            if component == None and (kb == None or ku == None or cooperativity == None):
                raise ValueError("Must pass in a Component or values for kb, ku, and coopertivity.")
            binder_params[binder] = {"kb":kb,"ku":ku,"cooperativity":coop_val}
        #out_rxns = []
        rxndict = {}
        coop_dict = {a.name:binder_params[a]["cooperativity"] for a in binder_params}
        for i in range(1, len(binders)+1):
            for combo in it.combinations(binders,i):
                #come up with all combinations of binders
                
                product = self.make_cooperative_complex(combo,bindee,coop_dict)
                #this is the complex which becomes the product
                for binder in combo:
                    #then in each case, one binder can be added to make
                    #this combination
                    reactant = tuple(set(combo)-set([binder]))
                    rxn_prototype = (binder,reactant)
                    
                    #this part makes a describer of the reaction; which reactants are combining?
                    #print(rxn_prototype)
                    #print(rxndict)
                    if(rxn_prototype in rxndict):
                        #if we already did this reaction then forget about it
                        continue
                    else:
                        reactant_complex = self.make_cooperative_complex(reactant,bindee,coop_dict)
                        reaction = Reaction(inputs=[binder, reactant_complex], outputs=[product],
                        input_coefs=[binder_params[binder]["cooperativity"], 1], output_coefs=[1], \
                                         k=binder_params[binder]["kb"],k_rev=binder_params[binder]["ku"])
                        rxndict[rxn_prototype]=reaction
        return [rxndict[a] for a in rxndict]
In [4]:
class CombinatorialPromoter(Promoter):
    def __init__(self, name, regulators, leak = True, assembly = None,
                 transcript = None, length = 0, mechanisms = {},
                 parameters = {},tx_capable_list = None,cooperativity = None, **keywords):
        """
        tx_capable_list = [[1,2,3],[0,2,3]] this means having regulator 1, 2, and 3 will transcribe
                                           but also 3, 2, and 0.
                                           #TODO make it force sorted list
        """
        if not isinstance(regulators, list):
            #you could give one string as a regulator
            regulators = [regulators]
        self.cooperativity = cooperativity
        self.regulators = []
        for regulator in regulators:
            if(isinstance(regulator,str)):
                self.regulators += [self.set_species(regulator, material_type = "protein")]
                #if it's a string then assume it's a protein
            elif(isinstance(regulator,Species)):
                #if it's already a species then add it wholesale
                self.regulators += [regulator]
        #after we've sanitized the inputs, then sort
        self.regulators = sorted(self.regulators)
        #now let's work out the tx_capable_list
        if(tx_capable_list == None):
            #if nothing is passed assume default
            self.tx_capable_list = [set([a.name for a in self.regulators])]
        elif(type(tx_capable_list)==list):
            #if the user passed a list then the user knows what they want
            self.tx_capable_list = [set(a) for a in tx_capable_list]

        self.leak = leak

        self.default_mechanisms = {"binding": Combinatorial_Cooperative_Binding()}

        Promoter.__init__(self, name = name, assembly = assembly,
                          transcript = transcript, length = length,
                          mechanisms = mechanisms, parameters = parameters,
                          **keywords)
        self.complex_combinations = {}
        self.tx_capable_complexes = []
    def update_species(self):

        mech_tx = self.mechanisms["transcription"]
        mech_b = self.mechanisms['binding']
        #set the tx_capable_complexes to nothing because we havent updated species yet!
        self.tx_capable_complexes = []
        species = []
        self.complexes = []
        if self.leak is not False:
            species += mech_tx.update_species(dna = self.assembly.dna)

        bound_species = mech_b.update_species(self.regulators,self.assembly.dna,\
                        component = self,part_id = self.name,cooperativity=self.cooperativity)
        #above is all the species with DNA bound to regulators. Now, we need to extract only the ones which
        #are transcribable
        for bound_complex in bound_species:
            species_inside = []
            for regulator in self.regulators:
                if(regulator.name in bound_complex.name):
                    species_inside += [regulator.name]
            #print("species_inside is "+str(species_inside))
            #print("tx capable is " +str(self.tx_capable_list))
            #print(set(species_inside) in [set(a) for a in self.tx_capable_list])
            if(set(species_inside) in [set(a) for a in self.tx_capable_list]):

                tx_capable_species = mech_tx.update_species(dna = bound_complex, transcript = self.transcript, \
                                                                                    protein = self.assembly.protein)
                species +=tx_capable_species[1:]
                self.tx_capable_complexes +=[bound_complex]
        species+=bound_species
        return species

    def update_reactions(self):
        reactions = []
        mech_tx = self.mechanisms["transcription"]
        mech_b = self.mechanisms['binding']

        if self.leak != False:
            reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
                                                            transcript = self.transcript, protein = self.assembly.protein)

        reactions += mech_b.update_reactions(self.regulators,self.assembly.dna,component = self,\
                                                        part_id = self.name,cooperativity=self.cooperativity)
        if(self.tx_capable_complexes == None or self.tx_capable_complexes == []):
            #this could mean we haven't run update_species() yet
            species = self.update_species()
            if(self.tx_capable_complexes == []):
                #if it's still zero after running update_species then we could be in trouble
                warn("nothing can transcribe from combinatorial promoter {}".format(self.name))

        else:
            for specie in self.tx_capable_complexes:
                tx_partid = self.name
                for part in specie.species_set:
                    #construct the name of the promoter with regulators bound
                    if part.material_type == "dna":
                        #the DNA doesn't matter
                        pass
                    else:
                        #put in the regulators!
                        tx_partid += "_"+part.name
                if(tx_partid[0]=="_"):
                    tx_partid = tx_partid[1:]
                #if it's bound to RNAP then it transcribes, right?
                tx_partid = tx_partid+"_RNAP"
                #print("tx_partid is "+ str(tx_partid))
                reactions += mech_tx.update_reactions(dna = specie, component = self, part_id = tx_partid, \
                                            transcript = self.transcript, protein = self.assembly.protein)

        return reactions

Tetrathionate Sensing Circuit in Bioscrape

Sensory kinase, ttrS, is membrane bound and will phosphorylate when tetrathionate is present. In turn, it phosphorylates ttrR. A ttrR-phosphorylated dimer is needed to activate the pTTR promoter.

TTRFALL.png

Let's create a mass action model for the system.

\begin{align} \textrm{Polymerase binds to P1d: } P + P1d &\underset{k^u_1}{\overset{k^b_1}\rightleftharpoons} P1d : P \\ \textrm{Transcription of ttrS, ttrR: } P1d : P &\xrightarrow{k_{tx}} ttrS_T + ttrR_T + P1d + P\\ \textrm{Translation ttrS_T to ttrS: } ttrS_T + R &\underset{k^u_2}{\overset{k^b_2}\rightleftharpoons} ttrS_T: R \xrightarrow{k_{tl}} ttrS_T + R + ttrS \\ \textrm{Translation ttrR_T to ttrR: } ttrR_T + R &\underset{k^u_3}{\overset{k^b_3}\rightleftharpoons} ttrR_T: R \xrightarrow{k_{tl}} ttrR_T + R + ttrR \\ \textrm{ttrS phosphorylation: } ttrS + tt &\underset{k^u_4}{\overset{k^b_4}\rightleftharpoons} ttrS : P + tt\\ \textrm{ttrR binding to unphosphorylated ttrS: } ttrR + ttrS &\underset{k^u_5}{\overset{k^b_5}\rightleftharpoons} ttrR:ttrS\\ \textrm{ttrR binding to phosphorylated ttrS: } ttrR + ttrS:P &\underset{k^u_6}{\overset{k^b_6}\rightleftharpoons} ttrR:ttrS:P\\ \textrm{Phosphorylation of ttrR: } ttrR:ttrS:P &\underset{k^u_7}{\overset{k^b_7}\rightleftharpoons} ttrR:P + ttrS\\ \textrm{Dephosphorylation of ttrR: } ttrR:P &\xrightarrow{k_{dephos}} ttrR \\ \textrm{Dimerization of ttrR:P: } 2ttrR:P &\underset{k^u_8}{\overset{k^b_8}\rightleftharpoons} ttrR:P_{dimer} \\ \textrm{ttrR:P dimer binding to pttr promoter: } ttrR:P_{dimer} + pttr0 &\underset{k^u_9}{\overset{k^b_9}\rightleftharpoons} pttr1 \\ \textrm{Polymerase binds to pttr: } P + pttr1 &\underset{k^u_{10}}{\overset{k^b_{10}}\rightleftharpoons} pttr1 : P \\ \textrm{Transcription of GFP: } pttr1 : P &\xrightarrow{k_{tx}} pttr1 + P + GFP_T \\ \textrm{Translation GFP: } GFP_T + R &\underset{k^u_{11}}{\overset{k^b_{11}}\rightleftharpoons} GFP_T: R \xrightarrow{k_{tl}} GFP_T + R + GFP \\ ttrS_T &\xrightarrow{\delta} \emptyset\\ ttrS &\xrightarrow{\delta} \emptyset\\ ttrR_T &\xrightarrow{\delta} \emptyset\\ ttrR &\xrightarrow{\delta} \emptyset\\ GFP_T &\xrightarrow{\delta} \emptyset \\ GFP &\xrightarrow{\delta} \emptyset \end{align}
In [5]:
#Species: A list of all the species
species = ["P", "P1d", "P1d:P", "ttrS_T", "ttrR_T", "R", "ttrR_T:R", "ttrS_T:R", "tt", "ttrR", "ttrS",
           "ttrS:P", "ttrR:ttrS", "ttrR:ttrS:P", "ttrR:P", "ttrR:P_dimer", "pttr0", "pttr1", "pttr1:P", "GFP_T", "GFP_T:R", "GFP"]
In [6]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 50), ("k2u", 10), 
          ("k3b", 30), ("k3u", 10), 
          ("k4b", 80), ("k4u", 2),
          ("k5b", 50), ("k5u",1), 
          ("k6b", 50), ("k6u", 1), 
          ("k7b", 50), ("k7u", 1), 
          ("k8b", 80), ("k8u", 1), 
          ("k9b", 50), ("k9u", 0),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 1),     
          ("ktx", .5), ("ktl", 5), ("k_dephos", 0.1), ("delta", .5)]
In [7]:
#Reactions: A list of reactions [rxn1, rxn2...]. Each reaction is a tuple ([Input Species], [Output Species], "propensity_type", {propensity_parameters})

# P + P1d <-> P1d:P
rxn1b, rxn1u = (["P", "P1d"], ["P1d:P"], "massaction", {"k":"k1b"}), (["P1d:P"], ["P", "P1d"], "massaction", {"k":"k1u"})

# P1d:P -> ttrS_T + ttrR_T
rxntx1 = (["P1d:P"], ["ttrS_T", "ttrR_T"], "massaction", {"k":"ktx"})

# ttrS_T + R <-> ttrS_T:R
rxn2b, rxn2u = (["ttrS_T", "R"], ["ttrS_T:R"], "massaction", {"k":"k2b"}), (["ttrS_T:R"], ["ttrS_T", "R"], "massaction", {"k":"k2u"})

# ttrS_T:R -> ttrS_T + R + ttrS
rxntl1 = (["ttrS_T:R"], ["ttrS_T", "ttrS", "R"], "massaction", {"k":"ktl"})

# ttrR_T + R <-> ttrR_T:R
rxn3b, rxn3u = (["ttrR_T", "R"], ["ttrR_T:R"], "massaction", {"k":"k3b"}), (["ttrR_T:R"], ["ttrR_T", "R"], "massaction", {"k":"k3u"})

# ttrR_T:R -> ttrR_T + R + ttrR
rxntl2 = (["ttrR_T:R"], ["ttrR_T", "ttrR", "R"], "massaction", {"k":"ktl"})

# ttrS + tt <-> ttrS:P + tt
rxn4b, rxn4u = (["ttrS", "tt"], ["ttrS:P", "tt"], "massaction", {"k":"k4b"}), (["ttrS:P", "tt"], ["ttrS", "tt"], "massaction", {"k":"k4u"})

# ttrR + ttrS <-> ttrR:ttrS
rxn5b, rxn5u = (["ttrR", "ttrS"], ["ttrR:ttrS"], "massaction", {"k":"k5b"}), (["ttrR:ttrS"], ["ttrR", "ttrS"], "massaction", {"k":"k5u"})

# ttrR + ttrS <-> ttrR:ttrS:P
rxn6b, rxn6u = (["ttrS:P", "ttrR"], ["ttrR:ttrS:P"], "massaction", {"k":"k6b"}), (["ttrR:ttrS:P"], ["ttrS:P", "ttrR"], "massaction", {"k":"k6u"})

# ttrR:ttrS:P <-> ttrR:P + ttrS
rxn7b, rxn7u = (["ttrR:ttrS:P"], ["ttrR:P", "ttrS"], "massaction", {"k":"k7b"}), (["ttrR:P", "ttrS"], ["ttrR:ttrS:P"], "massaction", {"k":"k7u"})

# ttrR:P -> ttrR
rxn_dephos = (["ttrR:P"], ["ttrR"], "massaction", {"k":"k_dephos"})

# 2ttrR:P <-> ttrR:P_dimer
rxn8b, rxn8u = (["ttrR:P", "ttrR:P"], ["ttrR:P_dimer"], "massaction", {"k":"k8b"}), (["ttrR:P_dimer"], ["ttrR:P", "ttrR:P"], "massaction", {"k":"k8u"})

# ttrR:P_dimer + pttr0 <-> pttr1
rxn9b, rxn9u = (["ttrR:P_dimer", "pttr0"], ["pttr1"], "massaction", {"k":"k9b"}), (["pttr1"], ["ttrR:P_dimer", "pttr0"], "massaction", {"k":"k9u"})

# P + pttr1 <-> pttr1:P
rxn10b, rxn10u = (["P", "pttr1"], ["pttr1:P"], "massaction", {"k":"k10b"}), (["pttr1:P"], ["P", "pttr1"], "massaction", {"k":"k10u"})

# pttr1:P -> pttr1 + P + GFP_T
rxntx2 = (["pttr1:P"], ["pttr1", "P", "GFP_T"], "massaction", {"k":"ktx"})

# GFP_T + R <-> GFP_T:R
rxn11b, rxn11u = (["GFP_T", "R"], ["GFP_T:R"], "massaction", {"k":"k11b"}), (["GFP_T:R"], ["GFP_T", "R"], "massaction", {"k":"k11u"})

# GFP_T:R -> GFP_T + GFP + R
rxntl3 = (["GFP_T:R"], ["GFP_T", "GFP", "R"], "massaction", {"k":"ktl"})


# RNA deg
rxnd_ttrS_T = (["ttrS_T"], [], "massaction", {"k":"delta"})
rxnd_ttrR_T = (["ttrR_T"], [], "massaction", {"k":"delta"})
rxnd_GFP_T = (["GFP_T"], [], "massaction", {"k":"delta"})

# Protein Deg
rxnd_ttrS = (["ttrS"], [], "massaction", {"k":"delta"})
rxnd_ttrR = (["ttrR"], [], "massaction", {"k":"delta"})
rxnd_GFP = (["GFP"], [], "massaction", {"k":"delta"})

reactions = [rxn1b, rxn1u, rxn2b, rxn2u, rxn3b, rxn3u, rxn4b, rxn4u,
             rxn5b, rxn5u, rxn6b, rxn6u, rxn7b, rxn7u, rxn8b, rxn8u,
             rxn9b, rxn9u, rxn10b, rxn10u, rxn11b, rxn11u,
             rxntx1, rxntx2,
             rxntl1, rxntl2, rxntl3,
             rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
             rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [8]:
#An initial condition for each species (uninitialized species default to 0)
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":25,
    "R":100,
}
In [9]:
#Instantiate the Model [The only object must of us will care about]
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
timepoints = np.linspace(0, 600, 100)

Results = py_simulate_model(timepoints, Model = M, stochastic = False) #py_simulate_model takes care of all other objects for you
Results.head()
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:2: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...
Out[9]:
P P1d P1d:P ttrS_T R ttrS_T:R ttrR_T ttrR_T:R ttrS tt ... ttrR:ttrS:P ttrR:P ttrR:P_dimer pttr0 pttr1 pttr1:P GFP_T GFP_T:R GFP time
0 25.000000 5.000000e+00 0.000000 0.000000 100.000000 0.000000 0.000000 0.000000 0.000000 200.0 ... 0.000000 0.000000 0.000000 5.000000e+00 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
1 15.010135 1.611876e-04 0.241864 0.018063 78.237244 4.703195 0.029872 4.667228 2.419410 200.0 ... 0.506453 0.843235 45.286372 -1.004820e-10 0.009973 4.990027 0.098135 12.392334 79.046695 6.060606
2 15.009981 7.801303e-06 0.011706 0.022982 63.516113 4.866145 0.037675 4.786917 5.689629 200.0 ... 0.622999 1.271489 117.406207 -1.967274e-16 0.009973 4.990027 0.257132 26.830824 221.174533 12.121212
3 15.009974 3.775699e-07 0.000567 0.028791 49.938064 4.793574 0.046662 4.662080 8.725189 200.0 ... 0.742075 1.583019 188.852833 -1.463558e-27 0.009973 4.990027 0.492288 40.606282 361.094857 18.181818
4 15.009973 1.827384e-08 0.000027 0.037507 37.469666 4.685892 0.059877 4.489338 11.479016 200.0 ... 0.869749 1.834576 258.058291 -3.773764e-25 0.009973 4.990027 0.859681 53.355103 492.248033 24.242424

5 rows × 23 columns

Matplotlib plotting of Bioscrape results

In [10]:
#Bioscrape returns a Pandas Dataframe by default, so plotting is easy!

plt.figure(figsize = (20, 10))
plt.subplot(221)
plt.title("Tetrathionate Regulators")
plt.plot(timepoints, Results["ttrR"], color = "pink", label = "ttrR")
plt.plot(timepoints, Results["ttrR_T"] + Results["ttrR_T:R"],":", color = "pink",  label = "ttrR_T")
plt.plot(timepoints, Results["ttrS"], color = "blue", label = "ttrS")
plt.plot(timepoints, Results["ttrS_T"] + Results["ttrS_T:R"],  ":", color = "blue",label = "ttrS_T")
plt.plot(timepoints, Results["tt"], color = "red", label = "tt")
plt.legend()

plt.subplot(222)
plt.title("Phosphorylation Pathway")
plt.plot(timepoints, Results["ttrR:ttrS"], label = "ttrR:ttrS")
plt.plot(timepoints, Results["ttrR:ttrS:P"], label = "ttrR:ttrS:P")
plt.plot(timepoints, Results["ttrR:P"], label = "ttrR:P")
plt.plot(timepoints, Results["ttrR:P_dimer"], label = "ttrR:P dimer")
plt.legend()

plt.subplot(223)
plt.title("Promoters")
plt.plot(timepoints, Results["P1d:P"], color = 'blue', label = "P1d:P")
plt.plot(timepoints, Results["P1d"], ":", color = 'blue', label = "P1d")
plt.plot(timepoints, Results["pttr1"], color = "cyan", label = "pttr1")
plt.plot(timepoints, Results["pttr0"], ":", color = "cyan", label = "pttr0")
plt.legend()

plt.subplot(224)
plt.title("Readout")
plt.plot(timepoints, Results["GFP_T"] + Results["GFP_T:R"], ":", color = 'green', label = "GFP mRNA")
plt.plot(timepoints, Results["GFP"], color = 'green', label = "GFP")

plt.legend()
Out[10]:
<matplotlib.legend.Legend at 0x1a2a6b6e50>

Holoviews plotting of Bioscrape results

Sort out Plotting

In [11]:
def plot_holoviews_ttr_bioscrape(results_melted):
    plot_group = []

    for row in results_melted.iterrows():
        sp = row[1]['species']

        if sp in ["P1d", "P1d:P", "pttr0", "pttr1", "pttr1:P"]:
            plot_group.append('promoters')

        elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T:R"]:
            plot_group.append('mRNA')

        elif sp in ["ttrS:P", "ttrR:ttrS", "ttrR:ttrS:P", "ttrR:P", "ttrR:P_dimer"]:
            plot_group.append('complexes')

        elif sp in ["ttrR", "ttrS", "GFP"]:
            plot_group.append('proteins')

        elif sp in ["P",   "R", "tt"]:
            plot_group.append('machinery and signals')
        else:
            print(sp)

    results_melted['plot_group'] = plot_group
    

    plot_list = []
    for group in results_melted['plot_group'].unique():
        plot_list.append(hv.Points(
            data=results_melted.loc[results_melted['plot_group'] == group],
            kdims=['time', 'value'],
            vdims = ['species']
        ).groupby(
            [ 'species']
        ).opts(
            tools=['hover'],
            legend_position = 'right',
            title = f'{group}',
            width = 500,
            #shared_axes=False
        ).overlay('species'))
        
    return plot_list

Now, proceed with plotting

In [12]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[12]:
time species value
0 0.000000 P 25.000000
1 6.060606 P 15.010135
2 12.121212 P 15.009981
3 18.181818 P 15.009974
4 24.242424 P 15.009973
In [13]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(5).opts(shared_axes=False)
Out[13]:

Decrease ttrS phosphorylation

Now, I am curious how the time course of the complexes changes if we decrease the sensory kinase's rate of phosphorylating. Perhaps the cell is low on phosphates due to other processes, and thus this rate is lower. We will thus just change the parameter associated with the forward reaction of ttrS phosphorylation.

In [14]:
M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)

# See what happens if our sensory molecule is slower
M.set_parameter("k4b",0.0001)

#Simulate the Model
timepoints = np.linspace(0, 600, 100)

Results = py_simulate_model(timepoints, Model = M, stochastic = False)
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  """Entry point for launching an IPython kernel.

We can use our plotting function we defined above:

In [15]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(5).opts(shared_axes=False)
Out[15]:

As expected, the entire response involves a delay now. Since phosphorylation is slower, we actually have no accumulation of ttrS:P. It is used immediately after it is created, and never builds up. We do have a lot of ttrR:ttrS compelxes, which are "stuck"

Test if Dose-Response Curve agrees with experimental data

In [16]:
M.set_parameter("delta", .5)

GFP_max = []

#Different initial values of S
tts = np.logspace(-3, 3, 100)

for tt0 in tts:
    x0["tt"] = tt0 #Change my initial condition dictionary
    M.set_species(x0)
    Results = py_simulate_model(timepoints, Model = M)
    GFP_max.append(np.max(Results["GFP"]))
In [17]:
LM_19 = pd.read_csv('LM19_data.csv')
LM_19
Out[17]:
Time variable OD Sample Concentration (mM) GFP (AU) GFP/OD700
0 23.95 A11 0.400 LM19 (high) 0.00 6688 16720
1 23.95 B11 0.373 LM19 (high) 0.10 38240 102520
2 23.95 C11 0.395 LM19 (high) 0.25 57234 144896
3 23.95 D11 0.393 LM19 (high) 0.50 71020 180712
4 23.95 E11 0.382 LM19 (high) 0.75 74950 196204
5 23.95 F11 0.379 LM19 (high) 1.00 75058 198042
6 23.95 G11 0.341 LM19 (high) 5.00 75993 222853
7 23.95 H11 0.318 LM19 (high) 10.00 73321 230569

Due to the unknown calibration factor between fluorescence and concentration, I have applied transformations to get the data in what seems to qualitatively match this curve. This is preliminary, and I hope to find a better way to show this agreement between data and simulation.

In [18]:
p = bokeh.plotting.figure(
    height=300,
    width=550,
    x_axis_label='[tt]',
    y_axis_label='Max GFP',
    x_axis_type = 'log',
    title = 'Comparison of Experimental Sigmoid with Simulated'
)

p.line(x=tts,y=GFP_max,)

# Try to transform data in some way that we can compare to simulation
p.circle(x = LM_19['Concentration (mM)'][1:]*4, y = LM_19['GFP/OD700'][1:]/150-600)

bokeh.io.show(p)

Changing Parameters to get ttr impact

Create Parameter Tabling Functions

In [19]:
rx_list = [
'Polymerase Bind to P1d',
'Ribosome binding to ttrS mRNA',
'Ribosome binding to ttrR mRNA',
'ttrS phosphorylation',
'ttR binds to unphosphorylated ttrS',
'ttR binds to phosphorylated ttrS',
'phosphorylation of ttrR',
'dimerization of ttrR~P',
'Dimer bindgint to pTTR',
'Polymerase binding to active pTTR',
'Ribosome binding to GFP mRNA',
    
]

extra_rx_list = ['transcription', 'translation', 'phosphorylation', 'degradation']
In [20]:
def print_rxn_df(parameters, x0):
    params_list_u = [p[0] for p in parameters if p[0].find('u') != -1]
    params_list_b = [p[0] for p in parameters if p[0].find('b') != -1]

    values_u = [p[1] for p in parameters if p[0].find('u') != -1]
    values_b = [p[1] for p in parameters if p[0].find('b') != -1]

    params_list_extra = [p[0] for p in parameters if p[0].find('b') == -1 and p[0].find('u') == -1]
    values_extra = [p[1] for p in parameters if p[0].find('b') == -1 and p[0].find('u') == -1]

    key_list = []
    val_list = []
    for key, value in x0.items() :
        key_list.append(key)
        val_list.append(value)
    
    params = pd.DataFrame()
    params['reaction'] = rx_list
    params['ku'] = values_u
    params['kb'] = values_b

    extra_params = pd.DataFrame()
    extra_params[''] = ['|']*len(values_extra)
    extra_params['extra reaction'] = extra_rx_list
    extra_params['k'] = values_extra

    inits = pd.DataFrame()
    inits[''] = ['|']*len(key_list)
    inits['initial'] = key_list
    inits['val'] = val_list

    return pd.concat([params, extra_params, inits],  axis=1).fillna("")

Plotting Function

In [21]:
def plot_ttr_curves(parameters, x0):
    experimental_ttr_concs = np.sort(LM_19['Concentration (mM)'][1:].values)
    r=list(bokeh.palettes.viridis(10))
    r.reverse()
    r = r[3:]

    M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
    timepoints = np.linspace(0, 1400, 1400)

    p = bokeh.plotting.figure(width = 500, height= 400, x_axis_label='Time (hr)', y_axis_label = 'GFP' )
    i = 0

    for tt_conc in tqdm.tqdm(experimental_ttr_concs):
        x0["tt"] = tt_conc #Change my initial condition dictionary
        M.set_species(x0)
        Results = py_simulate_model(timepoints, Model = M, stochastic = False)
        p.line(timepoints/60, Results['GFP'], color = r[i], legend_label = f'{tt_conc}', line_width = 2)
        i+=1

    p.legend.location = 'top_left'
    p.legend.title = 'Tetrathionate mM'

    return p

Inverse Tetrathionate Response

This set of parameter values gives higher activation when low tetrathionate as opposed to high tetrathionate.

In [22]:
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":250,
    "R":20,
}
In [23]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 0.004), ("k2u", 0.1), 
          ("k3b", 0.003), ("k3u", 5),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 0.0005), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 5), ("k5u",6), 
          ("k6b", 5), ("k6u", 10), 
          ("k7b", 1), ("k7u", 100), 
          ("k8b", 80), ("k8u", 5), 
          ("k9b", 400), ("k9u", 3),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 1),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 0.07), ("delta", .7)]
In [24]:
print_rxn_df(parameters, x0)
Out[24]:
reaction ku kb extra reaction k initial val
0 Polymerase Bind to P1d 1.0 100.0000 | transcription 0.1 | tt 200
1 Ribosome binding to ttrS mRNA 0.1 0.0040 | translation 1 | P1d 5
2 Ribosome binding to ttrR mRNA 5.0 0.0030 | phosphorylation 0.07 | pttr0 5
3 ttrS phosphorylation 0.5 0.0005 | degradation 0.7 | P 250
4 ttR binds to unphosphorylated ttrS 6.0 5.0000 | R 20
5 ttR binds to phosphorylated ttrS 10.0 5.0000
6 phosphorylation of ttrR 100.0 1.0000
7 dimerization of ttrR~P 5.0 80.0000
8 Dimer bindgint to pTTR 3.0 400.0000
9 Polymerase binding to active pTTR 1.0 50.0000
10 Ribosome binding to GFP mRNA 1.0 10.0000
In [25]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  import sys
100%|██████████| 7/7 [00:00<00:00, 96.09it/s]

Curves start together, then spread apart

In experiments, I found the different induction curves start close together and then fan out towards the end. I attempted find parameters that matched this behavior.

In [26]:
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":300,
    "R":10,
}
In [27]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 0.01), ("k2u", 100), 
          ("k3b", 0.3), ("k3u", 500),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 5), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 5), ("k5u",6), 
          ("k6b", 0.05), ("k6u", 1), 
          ("k7b", 0.0001), ("k7u", 1), 
          ("k8b", 8), ("k8u", 5), 
          ("k9b", 4), ("k9u", 3),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 10),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .7)]
In [28]:
print_rxn_df(parameters, x0)
Out[28]:
reaction ku kb extra reaction k initial val
0 Polymerase Bind to P1d 1.0 100.0000 | transcription 0.1 | tt 200
1 Ribosome binding to ttrS mRNA 100.0 0.0100 | translation 1 | P1d 5
2 Ribosome binding to ttrR mRNA 500.0 0.3000 | phosphorylation 7 | pttr0 5
3 ttrS phosphorylation 0.5 5.0000 | degradation 0.7 | P 300
4 ttR binds to unphosphorylated ttrS 6.0 5.0000 | R 10
5 ttR binds to phosphorylated ttrS 1.0 0.0500
6 phosphorylation of ttrR 1.0 0.0001
7 dimerization of ttrR~P 5.0 8.0000
8 Dimer bindgint to pTTR 3.0 4.0000
9 Polymerase binding to active pTTR 1.0 50.0000
10 Ribosome binding to GFP mRNA 10.0 10.0000
In [29]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  import sys
100%|██████████| 7/7 [00:00<00:00, 81.55it/s]

Curves start apart, then come together

Now, here is the case where induction curves start together and spread apart.

In [30]:
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":300,
    "R":25,
}
In [31]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 0.3), ("k2u", 10), 
          ("k3b", 0.03), ("k3u", 50),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 5), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 5), ("k5u",6), 
          ("k6b", 0.05), ("k6u", 1), 
          ("k7b", 0.0001), ("k7u", 1), 
          ("k8b", 8), ("k8u", 5), 
          ("k9b", 4), ("k9u", 3),
          ("k10b", 500), ("k10u",1), 
          ("k11b", 100), ("k11u", 100),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 70), ("delta", .7)]
In [32]:
print_rxn_df(parameters, x0)
Out[32]:
reaction ku kb extra reaction k initial val
0 Polymerase Bind to P1d 1.0 100.0000 | transcription 0.1 | tt 200
1 Ribosome binding to ttrS mRNA 10.0 0.3000 | translation 1 | P1d 5
2 Ribosome binding to ttrR mRNA 50.0 0.0300 | phosphorylation 70 | pttr0 5
3 ttrS phosphorylation 0.5 5.0000 | degradation 0.7 | P 300
4 ttR binds to unphosphorylated ttrS 6.0 5.0000 | R 25
5 ttR binds to phosphorylated ttrS 1.0 0.0500
6 phosphorylation of ttrR 1.0 0.0001
7 dimerization of ttrR~P 5.0 8.0000
8 Dimer bindgint to pTTR 3.0 4.0000
9 Polymerase binding to active pTTR 1.0 500.0000
10 Ribosome binding to GFP mRNA 100.0 100.0000
In [33]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  import sys
100%|██████████| 7/7 [00:00<00:00, 74.74it/s]

Amplified inverse response curve

Compared to the previous linear inverse response, the low induction curves grow in a concave curve in this example.

In [34]:
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":300,
    "R":1000,
}
In [35]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 0.3), ("k2u", 10), 
          ("k3b", 0.003), ("k3u", 50),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 5), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 500), ("k5u",6), 
          ("k6b", 0.05), ("k6u", 1), 
          ("k7b", 0.1), ("k7u", 1), 
          ("k8b", 8), ("k8u", 5), 
          ("k9b", 4), ("k9u", 3),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 10),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .7)]
In [36]:
print_rxn_df(parameters, x0)
Out[36]:
reaction ku kb extra reaction k initial val
0 Polymerase Bind to P1d 1.0 100.000 | transcription 0.1 | tt 200
1 Ribosome binding to ttrS mRNA 10.0 0.300 | translation 1 | P1d 5
2 Ribosome binding to ttrR mRNA 50.0 0.003 | phosphorylation 7 | pttr0 5
3 ttrS phosphorylation 0.5 5.000 | degradation 0.7 | P 300
4 ttR binds to unphosphorylated ttrS 6.0 500.000 | R 1000
5 ttR binds to phosphorylated ttrS 1.0 0.050
6 phosphorylation of ttrR 1.0 0.100
7 dimerization of ttrR~P 5.0 8.000
8 Dimer bindgint to pTTR 3.0 4.000
9 Polymerase binding to active pTTR 1.0 50.000
10 Ribosome binding to GFP mRNA 10.0 10.000
In [37]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  import sys
100%|██████████| 7/7 [00:00<00:00, 48.91it/s]

Close match to experiments

Now, we choose parameters that give decent quantitative agreement with the experimental data

In [38]:
x0 = {
    "tt":200,
    "P1d":5,
    "pttr0":5,
    "P":250,
    "R":20,
}
In [39]:
#Parameters: List of tuples [("paramter name" (string), value (float))]
parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 3), ("k2u", 10), 
          ("k3b", 0.00003), ("k3u", 5),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 5), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 5), ("k5u",6), 
          ("k6b", 5), ("k6u", 1), 
          ("k7b", 0.01), ("k7u", 1), 
          ("k8b", 80), ("k8u", 5), 
          ("k9b", 0.004), ("k9u", 3),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 1),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .9)]
In [40]:
print_rxn_df(parameters, x0)
Out[40]:
reaction ku kb extra reaction k initial val
0 Polymerase Bind to P1d 1.0 100.00000 | transcription 0.1 | tt 200
1 Ribosome binding to ttrS mRNA 10.0 3.00000 | translation 1 | P1d 5
2 Ribosome binding to ttrR mRNA 5.0 0.00003 | phosphorylation 7 | pttr0 5
3 ttrS phosphorylation 0.5 5.00000 | degradation 0.9 | P 250
4 ttR binds to unphosphorylated ttrS 6.0 5.00000 | R 20
5 ttR binds to phosphorylated ttrS 1.0 5.00000
6 phosphorylation of ttrR 1.0 0.01000
7 dimerization of ttrR~P 5.0 80.00000
8 Dimer bindgint to pTTR 3.0 0.00400
9 Polymerase binding to active pTTR 1.0 50.00000
10 Ribosome binding to GFP mRNA 1.0 10.00000
In [41]:
bokeh.io.show(plot_ttr_curves(parameters, x0))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 
  import sys
100%|██████████| 7/7 [00:00<00:00, 29.71it/s]

RBS Tuning

We take the above well-matching plot and perform RBS tuning, like we did experimentally with the ARL pool. We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively.

First, we will define a new plotting function that can take in different parameter values and plot on the same axis

In [42]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u):
    experimental_ttr_concs = np.sort(LM_19['Concentration (mM)'][1:].values)
    r=list(bokeh.palettes.viridis(10))
    r.reverse()
    r = r[3:]
    
    #Parameters: List of tuples [("paramter name" (string), value (float))]
    parameters = [("k1b", 100), ("k1u",1), 
          ("k2b", 3), ("k2u", k2u), 
          ("k3b", 0.00003), ("k3u", k3u),
          # Slow reverse reaction on ttrS phosphorylation
          ("k4b", 5), ("k4u", 0.5),
         # ttrR should quickly bind to unphos
          ("k5b", 5), ("k5u",6), 
          ("k6b", 5), ("k6u", 1), 
          ("k7b", 0.01), ("k7u", 1), 
          ("k8b", 80), ("k8u", 5), 
          ("k9b", 0.004), ("k9u", 3),
          ("k10b", 50), ("k10u",1), 
          ("k11b", 10), ("k11u", 1),     
          ("ktx", .1), ("ktl", 1), ("k_dephos", 7), ("delta", .9)]


    M = Model(species = [], reactions = reactions, parameters = parameters, initial_condition_dict = x0)
    timepoints = np.linspace(0, 1400, 1400)

    p = bokeh.plotting.figure(width = 500, height= 400, x_axis_label='Time (hr)', y_axis_label = 'GFP', title = f'ttrS (k2u): {k2u}, ttrR (k3u): {k3u}' )
    i = 0

    for tt_conc in experimental_ttr_concs:
        x0["tt"] = tt_conc #Change my initial condition dictionary
        M.set_species(x0)
        Results = py_simulate_model(timepoints, Model = M, stochastic = False)
        p.line(timepoints/60, Results['GFP'], color = r[i], legend_label = f'{tt_conc}', line_width = 2)
        i+=1

    p.legend.location = 'top_left'
    p.legend.title = 'Tetrathionate mM'
    p.y_range = Range1d(0,1.9*10**(-8))

    return p
In [43]:
# Define higher off rates
k2u = [20, 12, 7]
k3u = [15, 8, 5]

plot_list = []

for i in range(3):
    plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i]))
    
bokeh.io.show(gridplot([plot_list]))
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/ipykernel_launcher.py:24: UserWarning: The following species are uninitialized and their value has been defaulted to 0: P1d:P, ttrS_T, ttrS_T:R, ttrR_T, ttrR_T:R, ttrS, ttrS:P, ttrR, ttrR:ttrS, ttrR:ttrS:P, ttrR:P, ttrR:P_dimer, pttr1, pttr1:P, GFP_T, GFP_T:R, GFP, 

Great, this looks similar to the RBS tuning experiments we did.

Saving and Loading in as SBML

I did this as a check that SBML io ran properly.

In [44]:
# #Write the model M to a file
# #Note that bioscrape saves the current state of the model and will not auto-save any changes to it
# file_name = 'SignalTxTl.xml'
# M.write_bioscrape_xml(file_name) 

# M_loaded = Model(file_name)

# #M_SBML = import_sbml('sbml_file.xml') #This code won't work without first creating an SBML File

# #And the simulation will look the same as before
# Results = py_simulate_model(timepoints, Model = M_loaded, stochastic = False)

# plt.plot(timepoints, Results["S"]+Results["F1"], label = "S det", color = 'blue')
# plt.plot(timepoints, Results["T"]+Results["T:R"], label = "T det", color = 'cyan')
# plt.plot(timepoints, Results["X"], label = "X det", color = 'purple')
In [45]:
# M.write_sbml_model("TTR_only.xml")
In [46]:
M_ttr = import_sbml('TTR_only.xml')

#Simulate Deterministically and Stochastically
timepoints = np.linspace(0,700,10000)
result_det = py_simulate_model(timepoints, Model = M_ttr)
result_stoch = py_simulate_model(timepoints, Model = M_ttr, stochastic = True)
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/bioscrape/sbmlutil.py:208: UserWarning: Compartments, UnitDefintions, Events, and some other SBML model components are not recognized by bioscrape. Refer to the bioscrape wiki for more information.
  warnings.warn('Compartments, UnitDefintions, Events, and some other SBML model components are not recognized by bioscrape. ' +
In [47]:
import matplotlib as mpl
color_list = ['r', 'k', 'b','g','y','m','c']
mpl.rc('axes', prop_cycle=(mpl.cycler('color', color_list) ))
mpl.rc('xtick', labelsize=16) 
mpl.rc('ytick', labelsize=16)
plt.figure(figsize = (10, 4))
Out[47]:
<Figure size 720x288 with 0 Axes>
<Figure size 720x288 with 0 Axes>
In [48]:
#Plot Results
for i in range(len(M_ttr.get_species_list())):
    if i < len(color_list):
        s = M_ttr.get_species_list()[i]
        plt.plot(timepoints, result_det[s], color = color_list[i], label = "Deterministic "+s)
        plt.plot(timepoints, result_stoch[s], ":", color = color_list[i], label = "Stochastic "+s)

plt.title('TTR Model')
plt.xlabel('Time', FontSize = 16)
plt.ylabel('Amount', FontSize = 16)
plt.legend()
plt.show()

HRP AND Circuit

This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.

In [49]:
txtl = CRNLab("GFP")

txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')

#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")

dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")

#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()

hv.Points(
    data=results_melted.loc[results_melted['species'].str.startswith('protein')],
    kdims=['time', 'value'],
    vdims = ['species']
).groupby(
    ['species']
).opts(
    tools=['hover'],
    legend_position = 'right',
    width = 800,
    #shared_axes=False
).overlay()
Out[49]:

Vary hrpR and hrpS proteins

Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.

Get plotting in order

Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been

In [50]:
def simulate_heatmap(inducer_1_name,
                     inducer_1_values,
                     inducer_2_name,
                     inducer_2_values,
                     title,
                     x0=x0,
                     normalize=True,
                     max_to_1 = False,
                     CRN=CRN_extract_1,
                     timepoints=timepoints):
    
    GFP_max = pd.DataFrame()
    
    for conc_1 in inducer_1_values:
        x0[inducer_1_name] = conc_1 
        
        for conc_2 in inducer_2_values:
            x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
            
            Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
            GFP_max = GFP_max.append({inducer_1_name:conc_1,
                            inducer_2_name:conc_2,
                            'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)

    data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
                   columns = inducer_2_name,
                   values = 'GFP_max')
    
    fig, ax = plt.subplots()

    if normalize:
        data /= np.max(GFP_max['GFP_max'])
        
    if max_to_1:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.1, vmin = 0, vmax = 1)
    else:
        heat_map = sns.heatmap(data,
                         cmap = 'viridis',
                        linewidths=.2)

    heat_map.set_title(title)
    fig.tight_layout()
    return fig

Normal Functioning And gate

This non-leaky case shows that when we increase protein hrpR and hrpS, we see a continued increase in the amount of GFP produced.

In [51]:
hrpS_values = [0, 0.001, 0.01, 0.1, 0.5, 1, 10]
hrpR_values = [0, 0.0001, 0.001, 0.01, 0.1, 0.5, 1]

fig = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='hrpR, hrpS protein AND gate')
#fig.savefig("test.png")
In [52]:
pd.read_csv('parameters_and.txt', sep="\t", header=0)[49:]
Out[52]:
mechanism part_id param_name param_val comments
49 NaN phrpL_hrpR cooperativity 1.00 hrpR binds 1 time to the hrp promoter
50 NaN phrpL_hrpS cooperativity 1.00 hrpS binds 1 time to the hrp promoter
51 NaN phrpL_hrpR ku 10.00 Unbinding of hrpR from hrp promoter
52 NaN phrpL_hrpS ku 10.00 Unbinding of hrpS from hrp promoter
53 NaN phrpL_hrpR_hrpS_RNAP ktx 0.17 RNAP polymerase tx rate once it is bound to th...
54 NaN phrpL_hrpR_hrpS_RNAP ku 0.10 RNAP polymerase off rate once it is bound to t...
55 NaN phrpL_hrpR_RNAP ktx 0.15 Lower transcription rate if 1 activator
56 NaN phrpL_hrpR_RNAP ku 500.00 high unbinding rate when 1 activator
57 NaN phrpL_hrpS_RNAP ktx 0.15 Lower transcription rate if 1 activator
58 NaN phrpL_hrpS_RNAP ku 500.00 high unbinding rate when 1 activator

And Gate with leak in hrpR

Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.

In [53]:
pd.read_csv('parameters_and_leak.txt', sep="\t", header=0)[49:]
Out[53]:
mechanism part_id param_name param_val comments
49 NaN phrpL_hrpR cooperativity 1.00 hrpR binds 1 time to the hrp promoter
50 NaN phrpL_hrpS cooperativity 1.00 hrpS binds 1 time to the hrp promoter
51 NaN phrpL_hrpR ku 10.00 Unbinding of hrpR from hrp promoter
52 NaN phrpL_hrpS ku 10.00 Unbinding of hrpS from hrp promoter
53 NaN phrpL_hrpR_hrpS_RNAP ktx 0.17 RNAP polymerase tx rate once it is bound to th...
54 NaN phrpL_hrpR_hrpS_RNAP ku 0.10 RNAP polymerase off rate once it is bound to t...
55 NaN phrpL_hrpR_RNAP ktx 0.15 Lower transcription rate if 1 activator
56 NaN phrpL_hrpR_RNAP ku 5.00 high unbinding rate when 1 activator
57 NaN phrpL_hrpS_RNAP ktx 0.15 Lower transcription rate if 1 activator
58 NaN phrpL_hrpS_RNAP ku 500.00 high unbinding rate when 1 activator
In [54]:
txtl = CRNLab("GFP")

txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')


#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")


dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")

#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and_leak.txt")
#print("BioCRNpyler Representation:\n", repr(extract_1_TXTL))
CRN_extract_1 = extract_1_TXTL.compile_crn()
#CRN_model = extract_1_TXTL.get_model()
#print("\nCRN Representation:\n", repr(CRN_extract_1))

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()

hv.Points(
    data=results_melted.loc[results_melted['species'].str.startswith('protein')],
    kdims=['time', 'value'],
    vdims = ['species']
).groupby(
    ['species']
).opts(
    tools=['hover'],
    legend_position = 'right',
    width = 800,
    #shared_axes=False
).overlay()
Out[54]:
In [55]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]

fig = simulate_heatmap('protein_hrpR',hrpR_values,'protein_hrpS',hrpS_values,x0=x0,normalize=False,title='hrpR, hrpS protein AND gate', CRN=CRN_extract_1)
fig.savefig("leak.png")

And Gate with lots leak in hrpR becomes NAND

In [56]:
pd.read_csv('parameters_and_hrpR.txt', sep="\t", header=0)[49:]
Out[56]:
mechanism part_id param_name param_val comments
49 NaN phrpL_hrpR cooperativity 1.00 hrpR binds 1 time to the hrp promoter
50 NaN phrpL_hrpS cooperativity 1.00 hrpS binds 1 time to the hrp promoter
51 NaN phrpL_hrpR ku 10.00 Unbinding of hrpR from hrp promoter
52 NaN phrpL_hrpS ku 10.00 Unbinding of hrpS from hrp promoter
53 NaN phrpL_hrpR_hrpS_RNAP ktx 0.17 RNAP polymerase tx rate once it is bound to th...
54 NaN phrpL_hrpR_hrpS_RNAP ku 0.10 RNAP polymerase off rate once it is bound to t...
55 NaN phrpL_hrpR_RNAP ktx 150.00 Lower transcription rate if 1 activator
56 NaN phrpL_hrpR_RNAP ku 0.50 high unbinding rate when 1 activator
57 NaN phrpL_hrpS_RNAP ktx 0.15 Lower transcription rate if 1 activator
58 NaN phrpL_hrpS_RNAP ku 500.00 high unbinding rate when 1 activator
In [57]:
txtl = CRNLab("GFP")

txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = 'BasicExtract.tsv')


#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")


dna = DNAassembly("mydna",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")

#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "parameters_and_hrpR.txt")
#print("BioCRNpyler Representation:\n", repr(extract_1_TXTL))
CRN_extract_1 = extract_1_TXTL.compile_crn()
#CRN_model = extract_1_TXTL.get_model()
#print("\nCRN Representation:\n", repr(CRN_extract_1))

timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_mydna":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)

results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()

hv.Points(
    data=results_melted.loc[results_melted['species'].str.startswith('protein')],
    kdims=['time', 'value'],
    vdims = ['species']
).groupby(
    ['species']
).opts(
    tools=['hover'],
    legend_position = 'right',
    width = 800,
    #shared_axes=False
).overlay()
Out[57]:
In [58]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1, 10.0]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1, 10.0]

fig = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=False,title='hrpR, hrpS protein AND gate', CRN=CRN_extract_1)
fig.savefig("nand.png")
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...

AND gate

LacI promoter class

In [59]:
class LacIPromoter2(Promoter):
    def __init__(self, name = "pLac",  IPTG = "IPTG",
                 lacI = "lacI", leak = True, assembly = None,
                 transcript = None, length = 0, mechanisms = {},
                 parameters = {}, **keywords):

        Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
        #RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)

        # Transcription is already included in promoter class
        self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}

        self.leak = leak

        self.IPTG = self.set_species(IPTG, material_type = "inducer")
        self.lacI = self.set_species(lacI, material_type = "protein")

        self.plac_2_lacI = None
        self.IPTG_2_lacI = None

    def update_species(self):
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']
        species = []
        self.complexes = []

        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + pLac <--> 2xlacI:pLac
        """
        species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
        species += species_b

        #print('plac stuff \n',species_b)
        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
                self.plac_2_lacI = s

        if self.leak is not False:
            species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        # IPTG can bind to lacI to sequester it
        # give citation for 1 being enough
        """A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
           2lacI <--> 2xlacI
           2xlacI + IPTG <--> 2xlacI:IPTG
        """
        species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
        #print('iptg stuff \n', species_b)

        species += species_b

        for s in species_b:
            if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
                self.IPTG_2_lacI = s

        # Unbound transcription
        species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return species


    def update_reactions(self):
        print('r')
        mech_tx = self.mechanisms['transcription']
        mech_b = self.mechanisms['binding']

        if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
            self.update_species()

        if self.leak != False:
            # Will need transcription, plac_2x_lacI, ktx, ku, kb
            reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        reactions = []

        # Repression
        # Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
        reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
                                                                part_id = self.plac_2_lacI.name, cooperativity = 2)


        # LacI binding to IPTG
        # Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
        # .name gets rid of the complex
        # .type will give you complex, species, etc
        # .attributes gives you properties like phosphorylation state
        reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)

        '''
        Any bit of IPTG destroys the bound structure
        IPTG + plac_2xlacI --> IPTG:2xlacI + plac
        '''
        if isinstance(mech_b, Two_Step_Cooperative_Binding):
            kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
        else:
            kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)

        reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))

        # Unbound transcription plac, kb, ku, ktx
        reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
                                        transcript = self.transcript, protein = self.assembly.protein)

        return reactions

Sort out param printer

In [60]:
def param_printer2(x0, parameters):
    inits = pd.DataFrame(list(x0.items()),columns = ['param','value']) 
    

    df = pd.DataFrame(list(parameters.items()),columns = ['param','value']) 
    col_part_id = []
    mech = []
    for val in df['param']:
        if isinstance(val, tuple):
            col_part_id.append(val[0])
            mech.append(val[1])
        else:
            col_part_id.append(' ')
            mech.append(val)

    df['param'] = mech
    df['part'] = col_part_id
    df.reindex(['part','param','value'],axis=1)
    df[''] = ['|']*len(col_part_id)
    
    return pd.concat([df, inits],  axis=1).fillna("")

AND simulation

In [61]:
timepoints = np.linspace(0, 200, 1000)

x0 = {"dna_pttr_hrpR":5,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":10,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":5,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":10,
      "protein_Ribo":90.,
      "ligand_tt": 100,
      "inducer_IPTG": 10}
In [62]:
parameters = {
    "kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":2, 'cooperativity':2, #some default parameters
    
    ('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0, 
    ('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025,  #this is the transcription rate
    ('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('translation_mm', 'B0030', 'ku'): 10.0, #Unbinding
    ('translation_mm', 'B0030', 'ktl'): 1.5, #Translation Rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm','kb'):100, #default binding rate!
    
    #('ttrS-P',"kdephos"):1000, #auto dephosphorylation
    ('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
    
    #('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
    ('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
    
    #('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
    ('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
    ('ttrR-P',"kdephos"):1000, #auto dephosphorylation
    
    ('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
    ('phrpL_hrpS', 'cooperativity'): 1.0, 
    ('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025,  #this is the transcription rate
    ('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate

    ('translation_mm', 'B0030', 'ku'): 10.0, #Unbinding
    ('translation_mm', 'B0030', 'ktl'): 1.5, #Translation Rate

    ('transcription_mm',  'ktx'): .01,    #These are the parameters for transcription leak
    ('transcription_mm', 'ku'): 100,
    
    ('phrpL_hrpR','ku'):100,
    ('phrpL_hrpS','ku'):100,
    ('phrpL_hrpS_RNAP', 'ku'): 100.,#we'll say that leak has 10x the unbinding rate
    ('phrpL_hrpR_RNAP', 'ku'): 100,
    ('phrpL','ku'):100,
    
    ('binding', 'pLac_lacI','ktx'):0.001,
    ('plac_lacI','kb'):10,
    ('plac_lacI','ku'):0.001,
    
}


#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")

# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")

lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")

# # constit_P
# plac_hrpS = DNAassembly(name = "p_hrpS",promoter=RegulatedPromoter("P",["lacI"]),rbs="rbsH",protein="hrpS")


# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")

hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")

extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
                                                                     constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)

CRN_extract_1 = extract_1_TXTL.compile_crn()
r
In [63]:
param_printer2(x0, parameters)
Out[63]:
param value part param value
0 kb 100.000000 | dna_pttr_hrpR 5
1 ku 10.000000 | dna_constit_ttrS 1
2 ktx 0.050000 | dna_constit_ttrR 1
3 ktl 0.200000 | dna_constit_lacI 10
4 kdeg 2.000000 | dna_plac_hrpS 1
5 cooperativity 2.000000 | dna_hrpL_gfp 5
6 cooperativity 2.000000 pttr_phosphate_P_protein_ttrR | dna_constit_mScarlet 5
7 ktx 0.170250 pttr_phosphate_P_protein_ttrR | phosphate_P 100
8 ku 11.013216 pttr_phosphate_P_protein_ttrR | protein_RNAP 10
9 B0030 10.000000 translation_mm | protein_Ribo 90
10 B0030 1.500000 translation_mm | ligand_tt 100
11 ktx 0.010000 transcription_mm | inducer_IPTG 10
12 kb 100.000000 transcription_mm |
13 kdephos 1000.000000 ligand_tt_protein_ttrS-P |
14 kphos 10.000000 ligand_tt_protein_ttrS-P |
15 ktransfer 5.000000 ligand_tt_protein_ttrS-ttrR |
16 kdephos 1000.000000 ttrR-P |
17 cooperativity 1.000000 phrpL_hrpR |
18 cooperativity 1.000000 phrpL_hrpS |
19 ktx 0.170250 phrpL_hrpR_hrpS_RNAP |
20 ku 11.013216 phrpL_hrpR_hrpS_RNAP |
21 ku 100.000000 transcription_mm |
22 ku 100.000000 phrpL_hrpR |
23 ku 100.000000 phrpL_hrpS |
24 ku 100.000000 phrpL_hrpS_RNAP |
25 ku 100.000000 phrpL_hrpR_RNAP |
26 ku 100.000000 phrpL |
27 pLac_lacI 0.001000 binding |
28 kb 10.000000 plac_lacI |
29 ku 0.001000 plac_lacI |
In [64]:
#print("\nCRN Representation:\n", CRN_extract_1.pretty_print())
In [65]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)

fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
fig.savefig("and_first.png")
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...

Weird nand behavior again...

In [66]:
x0 = {"dna_pttr_hrpR":5,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":10,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":5,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":100,
      "protein_Ribo":90.,
      "ligand_tt": 100,
      "inducer_IPTG": 10}
In [67]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)

fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
#fig.savefig("and_first.png")
odeint failed with mxstep=500...odeint failed with mxstep=500...

Increasing RNAP, Ribo 10x from initial and

In [68]:
x0 = {"dna_pttr_hrpR":5,
      "dna_constit_ttrS":1,
      "dna_constit_ttrR":1,
      "dna_constit_lacI":10,
      "dna_plac_hrpS":1,
      "dna_hrpL_gfp":5,
      "dna_constit_mScarlet":5,
      "phosphate_P":100,
      "protein_RNAP":100,
      "protein_Ribo":900.,
      "ligand_tt": 100,
      "inducer_IPTG": 10}
In [69]:
param_printer2(x0, parameters)
Out[69]:
param value part param value
0 kb 100.000000 | dna_pttr_hrpR 5
1 ku 10.000000 | dna_constit_ttrS 1
2 ktx 0.050000 | dna_constit_ttrR 1
3 ktl 0.200000 | dna_constit_lacI 10
4 kdeg 2.000000 | dna_plac_hrpS 1
5 cooperativity 2.000000 | dna_hrpL_gfp 5
6 cooperativity 2.000000 pttr_phosphate_P_protein_ttrR | dna_constit_mScarlet 5
7 ktx 0.170250 pttr_phosphate_P_protein_ttrR | phosphate_P 100
8 ku 11.013216 pttr_phosphate_P_protein_ttrR | protein_RNAP 100
9 B0030 10.000000 translation_mm | protein_Ribo 900
10 B0030 1.500000 translation_mm | ligand_tt 100
11 ktx 0.010000 transcription_mm | inducer_IPTG 10
12 kb 100.000000 transcription_mm |
13 kdephos 1000.000000 ligand_tt_protein_ttrS-P |
14 kphos 10.000000 ligand_tt_protein_ttrS-P |
15 ktransfer 5.000000 ligand_tt_protein_ttrS-ttrR |
16 kdephos 1000.000000 ttrR-P |
17 cooperativity 1.000000 phrpL_hrpR |
18 cooperativity 1.000000 phrpL_hrpS |
19 ktx 0.170250 phrpL_hrpR_hrpS_RNAP |
20 ku 11.013216 phrpL_hrpR_hrpS_RNAP |
21 ku 100.000000 transcription_mm |
22 ku 100.000000 phrpL_hrpR |
23 ku 100.000000 phrpL_hrpS |
24 ku 100.000000 phrpL_hrpS_RNAP |
25 ku 100.000000 phrpL_hrpR_RNAP |
26 ku 100.000000 phrpL |
27 pLac_lacI 0.001000 binding |
28 kb 10.000000 plac_lacI |
29 ku 0.001000 plac_lacI |
In [70]:
IPTG_conc_list = np.logspace(-2, 4, 15).round(decimals=6)
tt_conc_list = np.logspace(-2, 4, 15).round(decimals=6)

fig = simulate_heatmap('ligand_tt',tt_conc_list,'inducer_IPTG',IPTG_conc_list,x0=x0,normalize=False,title='Compiled AND gate', CRN=CRN_extract_1)
#fig.savefig("resource.png")
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/chemical_reaction_network.py:870: UserWarning: The following species are uninitialized and their value has been defaulted to 0: inducer_IPTG, 
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...odeint failed with mxstep=500...

Computing Environment

In [71]:
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews
CPython 3.7.6
IPython 7.13.0

numpy 1.18.1
pandas 1.0.3
scipy 1.4.1
bioscrape 0.0.0
biocrnpyler 0.1
bokeh 2.0.1
jupyterlab 2.1.0
holoviews 1.13.2
/Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/watermark/watermark.py:157: UserWarning: Module biocrnpyler was already imported from /Users/lianamerk/anaconda3/envs/python3.7-bioscrape/lib/python3.7/site-packages/biocrnpyler-0.1-py3.7.egg/biocrnpyler/__init__.py, but /Users/lianamerk/git/BioCRNPyler is being added to sys.path
  ver = pkg_resources.get_distribution(p).version